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n^ ' We investigated numerically the dust growth driven by Brownian motion in a proto- 

' planetary disc around a solar-type young stellar object. This process is considered 

Cn ' as the first stage in the transformation of the initially micron-sized solid particles 

to a planetary system. In contrast to earlier studies the growth was investigated at 
the small particle number densities typical for the conditions in a proto-planetary 
disc. Under such circumstances, the mean particle distance exceeds the typical 
aggregate diameter by orders of magnitude, and a collision will be a very rare 
event. We derived a criterion which allows an efficient detection of candidates for 
Q^ ^ imminent collisions. The N-particle-method we used is based upon an adaptive 

time step scheme respecting the individual dynamical states of the aggregates. Its 
^sj ■ basic concept is to perform on average constant "length steps", instead of using 

r~^ ' constant time steps. The numerical cost of the algorithm scales with the particle 

/■^«, ' number better than N log N . In order to minimise the influence of the decreasing 



o\ 



number of particles within the simulation box, a new rescaling method is used 
throughout the aggregation process. Our numerical results indicate that at very 

r^ . low number densities, the growth process is influenced by spatial number density 

O ,■ fluctuations. 

o, 

H 1 Introduction 

^ ' Since the advent of high performance computers, scientists have been able to 

numerically investigate complex astronomical phenomena that have interested 
i^ astronomers and physicists for centuries. The first problems in numerical astro- 

I ^S , physics were hydrodynamical and radiative transfer problems. Then, starting 

C^ ' in the early 1980s, several workers began modeling-.the gravitational clustering 

of stars and galaxies efficiently using tree codestZl. Since the beginning of 
the 90' Particle codes are although used to investigate smaller-scale astrophysi- 
cal problems, such as the formation of the precursors of the planets (planetesimals). 

Nowadays it is a generally accepted view among the astrophysicists, that the 
genesis of a planetary system coincides with the formation of sun-like young stellar 
objects surrounded by a gaseous disc. The building bricks of the planetesimals 
are micron-sized solid particles (the so-called cosmic dust) embedded in the gas 
of the disc. The formation of a planetary system is a surprisingly fast process. 
From both, the astronomical observations of proto-planetary discs and geophysical 
studies, one can derive a formation time scale of a few million years. Up to the 
mid-90s the majority of the astrophysicists were convinced that the transition 
of the micron-sized dust_gtains to the kilometre-sized planetesimals were due to 
gravitational instabilitiesE3EI. So it was quite a surprise for the community, when 
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Weidenschillin^a and Cuzzi et al.H demonstrated that the inherent gas turbulence 
prevents dense dust layers from collapsing into a planetesimal. Another possible 
process for forming a planetesimal in a proto-planetary disc is the dust growth due 
to collisional sticking - the so-called coagulation. For particles to collide and stick, 
there must be a relative velocity component between the grains. In the onset of 
grain growth Brownian motion dominates other motionsEI. However, it can not be 
the source of an effective grain growth leading to planetesimals, because the growth 
time scales are much longer than the few million years. Consequently, there has 
to be a much more effective mechanism for driving the dust jCpagulation. Such a 
source of relative velocities could be the turbulent gas motiouEa if the typical time 
Tf a dust aggregate needs to couple with t|hO| gas stream, is larger than the typical 
life time of the smallest turbulent eddyElHl3. If one-considers realistic values for 
a proto-planetary disc around a solar-like protostarlil, then one finds that the Tf 
of the initial particles must increa|Sfi -bi^ five orders of magnitude in order to have 
grain growth driven by turbulencellj'EJ. Consequently, during the Brownian stage 
of dust growth the typical friction time of the dust aggregates t^it^j grow rapidly. 
However, Monte-Carlo type simulations of the cluster growthEjaO indicate that 
the growth tends to form fractal aggregates characterised by an almost constant 
coupling time Tf. In these simulations, particles are added to the growing dust 
grain according to some given rule. The disadvantage of such Monte-Carlo methods 
is that they require significant simplification of the astrophysical conditions. In 
addition, the Monte Carlo approach cannot address the inherently dynamic nature 
of the coagulation process. 

To overcome the disadvantages of Monte-Carlo astrophysical simulations, we 
developed a self-consistent N-particle code. In order to obtain a detailed under- 
standing of the coagulation in proto-planetary discs, we considered the infiuence of 
the aggregate structure in the coagulation process. The "classical" tool for a numer- 
ical study of this phenomenon is the molecular dynamic (MD). However, since the 
dust grains are performing a non-deterministic, diffusive motion, it is not straight- 
forward to apply the MD concept for modeling this process. Another difficulty 
with the MD approach is that we wish to study the coagulation at very low particle 
number densities. Or more precisely, the mean particle distance (x) = {3/4?™}^''^ 
(n - particle number density) exceeds the typical aggregate size by at least four 
orders of magnitude. Under such circumstances, a mutual collision of two particles 
is a rare event. In order to model the coagulation exactly, it is essential that no 
collisions are missed. To meet this requirement one needs a good estimate for the 
particle position after a certain period of time. In the following we demonstrate 
that, for sufficient small time steps, we can use the MD methods even for diffusing 
particles. 



2 Computational technique 
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2. 1 Simulation of diffusive trajectories 

Since the motion of growing dust grains is only affected by the gas-grain interaction 
(diffusion), and in some cases by the gravitational field, the particle dynamics be- 
tween successive collisions can be treated independently. The "equation of motion" 
for diffusive particles is given by the Langevin equation 

dtx — V 

dtv ^ -Tj^v + m^^Ffi. (f) 

Here, x and v are the Cartesian components of the position and velocity of the 
particle. Furthermore, the friction time r/ is the mean interval the particle needs 
to dissipate its relative kinetic energy with respect to the suspending gas. The 
influence of the gas-grain interaction is approximated by the stochastic force Ffi 
having the properties {Ffi) = and {F?i) ^ 2mkT/Tf (m - particle mass, T - 
kinetic gas temperature, k - Boltzmann constant). 

The diffusive path of the particle can be obtained by successive integration of 
(yJL^fter a time interval r — St/rf, the new particle position and velocity are given 
by@Q 

X = xo + T/(w + uo)tanh(T/2) + {2t^{v^)[t - 2 tanh (t/2)]}^^^ ^ (2) 

V = v„e-^ + {{v')[l - e-^^]Y^' ^2, (3) 

where ^i and ^2 are normalised Gaussian random numbers, i.e. {^i) — 0, (C1C2) — 0, 
and (^■^) = f (i = I72). As mentioned before, the MD concept requires knowledge 
of the particle position after a period of time. For a particle following a diffusive 
trajectory, it is impossible to obtain such an expression for an arbitrary interval dt, 
since only the probability of the particle being at position x after St can be derived. 
However, on time scales 5t<^Tf, the particle dynamics is still strongly correlated 
with its initial state. Performing a first order series expansion of Eq. (0) and (0) 
yields 

x = xn+vo5t + 0{T^'^) (4) 

v^vo[1-t] + ^21^(, + 0{t^'^). (5) 

Although Eq. (H) describes the motion of a free particle following a ballistic 
trajectory, the particle is moving stochastically for r ^ 1. Due to the "stochastic 
force" ^ in Eq. (0), the motion is not reversible in time. For this reason, stochastic 
path segments described by Eq. (0) and (ph are called pseudo ballistic. 

To obtain a test criterion for particle collisions, we make use of the fact that on 
time scales r ^ 1 the deviation of the actual trajectory from the ballistic path is 
small. From the higher terms of Eq. (0), one finds, that on average the dispersion 
of the particle position is 

((a; _ ^0 _ y^stf) ^ (Ax'') - l{v')Ty + 0{t^). (6) 

The dispersion in the direction of the ballistic trajectory (^ is negligible, while 
perpendicular to the trajectory, (Aa;' ) contributes significantly. Heuristically, this 
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result can be interpreted that the particle is moving towards its ballistic path (0) 
within a cone having the time-dependent opening radius (Aa;' )^''^ . Consequently, 




Figure 1. Example for two particles performing a pseudo-ballistic motion. A mutual collision can 
only happen with significant probability within the intersection of the cones. 



the prerequisite fpij an aggregate collision within t <C 1 is the intersection of their 
probability conestZl. If a particle pair fulfills this criterion, the trajectories of the 
two aggregates will be simulated with a high temporal resolution according to 
Eq. (y) and (|5|). the time step interval is obtained by the requirement that the 
faster particle move one diameter of the larger aggregate during that time. If the 
particles become closer to each other than the distance of the sum of their radii, 
then they will be combined into a new cluster. 

2.2 Numerical realisation: Calendar algorithm 

As mentioned before the coagulation of particles following a diffusive trajectory 
can be represented by MD methods if the applied time step 5t is much smaller 
than the friction time Tf . In the course of a time step St <l^ Tf, only particles closer 
than a critical separation are likely to collide. Only for these pairs the test of the 
collision criterion is required. Therefore, at each time step the nearest neighbour 
of the considered particle has to be determined. Experience shows that, due 
to the large velocity dispersion of the growing ensemble, algorithms with global 
time step schemes are inefficient for simulating the coagulation. Therefore we 
utilised a time step scheme respecting the dynamical state of the particles. The 
algorithm is similar to MD methods simulating hard sphere fluidal3. The basic 
idea is to perform average constant length steps, instead of constant time steps. 
To implement this idea, the simulation volume is divided into cubical cells of the 
same size. Next, the individual time step 5ti for a particle is defined by the interval 
Sti = t[—ti that the particle needs to move from its current cell to a neighbouring 
one, where ti and t[ are the instant of the last update of its coordinates (its "local 
time") and the instant of crossing the cell boundary, respectively. Obviously, 
within the interval (i^ . . . i^), a particle can only collide with particles also belonging 
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to its cell. It is essential that for the ensemble particle i leaving its cell first (ie 
t'j — mm{t'i . . -t'ff)), all possible collision partners have moved into its cell before 
ti. This leads to the following procedure: First, for each cell, build a list of the 
embedded particles. Next, for the particle with the smallest cell crossing time t^, 
the collision criterion is checked with respect to the other particles in its cell. If 
no collision is found, the particle position and velocity are updated according to 
Eq. (^ and (||). Finally, the cell list is updated and the procedure is repeated 
for the next particle with the momentarily smallest t'j (see Fig. 0). It should be 
noted that in contrast to "classical" MD, touching the cell boundary at t'^ does 
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Figure 2. A 2-dimcnsional example of three time steps according to the presented event-driven 
algorithm. The arrows show the predicted trajectories within the cells, while the italic numbers 
are the predicted instants t'^ of crossing the cell boundaries. For each time step the trajectory of 
the bold-printed particle is simulated. Before the first time step all particles have the local time 
ti = (Fig. a). Due to Stu = min{tj . . .^12) = 0.3, the first time step has to be conducted for 
particle 11. After this, the new local time for particle 11 is tii = 0.30 (Fig. b), while after the 
second time step (54i2 = min(i'j^ . . ■t'^2) = 0.35) the local time for particle 12 is ii2 = 0.35 (Fig. 
c). Note that for particle 12 periodic boundary conditions were used. 



not necessarily mean that the particle will instantaneously leave its embedding 
cell. Due to the stochastic nature of its motion, the particle can be "reflected" 
on the boundary. For a npre detailed consideration of such effects we refer to the 
description of our methodEil. 

The prerequisite for an efficient numerical implementation of the described pro- 
cedure is a fast search algorithm for the smallest cell crossing time t'^ (the so-called 
"event"). Therefore, the cell-crossing events are arranged in a tree-like structure 
(the so-called "calendar") as demonstrated in Fig. ||. 

2.3 Rescaling of the particle number 

Every N-particle method suffers from finite size effects due to the restricted number 
of simulated particles N . In order to minimise such effects, one chooses N as large 
as possible. However, for aggregation studies it is not sufficient just to start with a 
large N, since the particle number density decreases due to the coUisional sticking. 
In previous aggregation studies, the growth was simulated in a box of constant 
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Figure 3. Evolution of the event tree of the particle ensemble illustrated in Fig. bl The node 
colours correspond with the particle colours in Fig. H. In the upper part of the node the node 
index (identical with the particle index) as well as the predicted instants of cell crossing are given. 
The lower left number is the index of a linked later event, while the lower right number is the 
index of a linked earlier event. The bold arrows mark the search path to the earliest event, while 
the broken arrows identify the search path for sorting in a new event. 



volume. In that case, the total number N of simulated particles in the course of the 
simulation decreased, leading to strong finite size effects—Jn extreme cases, the box 
contained eventually only a single huge particle (eg ZiflnJ), which is an unphysical 
result. Such effects can be avoided if the number of simulated aggregates is kept 
constant on average. For maintaining an almost constant (N) we rescale the size 
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Figure 4. Two-dimensional example for the "refreshment" of the particle loss AN = N{t) — (N) 
by rescaling the box size. Fig. (a) shows the box before the rescaling. The hatched parts mark 
the areas copied to the new volume as demonstrated in Fig. (b) . On the average the rescaled box 
contains again (A') particles (Fig. (c)). 



L of the box according to X -> X' = L ^1 + AN{t)/{N} if the particle loss A A^ 
exceeds a given threshold. The new volume is refilled with copies of clusters of the 
original box (Fig. 0). Note, that the global observables of the simulation (mass 
density, cluster size spectrum) are not changed by this rescaling method. 
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3 Results 

In this section we give a "h r ipf overview of the results we obtained from simulations 
using scalar machinescSEj. We investigated the coagulation driven by Brownian 
motion for the typical conditions in a proto-planetary disc at 10 AU^ All simulation 
were started with N randomly distributed spheres of 1 micron diameter. The 
simplification of an initially mono-disperse size distribution is justified by the fact 
that the growth driven by Brownian motion tends to form aggregates of similar size. 
In order to study the influence of the initial dust number density no simulations 
with no between 10^ . . . 10^m~^ were performed. Most of the simulations were 
carried out with 10 000 particles. For each parameter set at least six runs with 
different seeds of the random generator were executed. 

3. 1 Structure of the growing aggregates 

As outlined in Sec. ^ the prerequisite for an effective grain growth driven by 
turbulence is a rapid increase of the aggregate friction time tj . Both, theoretical 
studiesQ and experimentso show, that for particles, which are small compared to 
the mean free path of the gas, t/ is proportional to their mass-surface ratio. This 
has the consequence, that the T/-evolution is determined by the structure of the 
growing aggregates. The increase of the particle mass m with respect to a typical 
radius i? is a characteristic scaling property of any growth regime and can be 
expressed by a fractal dimension Dm'- m{\R) — \^"^m{R). 

All published models for the dust coagulation assume that there is a unique 
value of the fractal dimension £>„ describing all aggregates of the dust ensemble. 
However, real growth processes produce a wide variety of different particle struc- 
tures. In (Fig. ^) the distribution of the fractal dimension of a simulated dust 
ensemble for fixed times is shown. We found that the mean fractal dimension {Dm) 
is about 1.8. This value is in good agreement with the results of Monte-Carlo simu- 
lations of the ballistic Cluster-Cluster- Aggregations (BCCA). In order to analyse 
the Tf evolution also the scaling of the projected aggregate area {a) with the typical 
radius has to be known. Our studies show (Fig. ^b) that the often used assumption 
of (7 ~ i?^ implies the wrong picture of the 7/ evolution. We found that a scales 
approximately as i?^'^, corresponding to a friction approximately proportional to 
the eighth root of the typical aggregate mass (Fig. ||c). Although we obtained in 
contrast to simplier models a size-dependent friction time our result still shows 
that growth solely driven by Brownian motion cannot form aggregates with friction 
times comparable to the lifetime of the smallest turbulent eddies within a million 
years. However, due to the broad distribution of the fractal dimension there are 
some aggregates growing as compact particles. These particles are of great im- 
portance because they could be possible seed grains for a runaway growth due to 
differential sedimentation. 



" We considered a disc of 80 AU with a mass accretion of M = lO~^M0yr~^ and an "alphap-i 
viscosity of o = 10~^. Using the one-dimensional hydrodynamical disc model of Ruden&PoUacko 
we determined the gas temperature Tg and pressure pg at 10 AU as about 100 K and 2 10^* Pa, 
respectively. 
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(a) scaling of the cluster mass: D (b) scaling of the projected area: D 
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Figure 5. Normalised distribution of the scaling exponents of the cluster mass (Fig. a), projected 
cluster area (Fig. b), and friction time (Fig. d) of all clusters consisting of more than 20 single 
spheres after 355 days of growth. The initial dust number density of the monomers was no = 
2-10^m~'^. Fig. c shows the mean size and size variation of the clusters contributing to the bin. 



3.2 Evolution of the mass spectrum 

In the previous section we discussed the friction time growth with the aggregate 
size. Equally important for the understanding of the formation of larger aggregates 
in proto-planetary discs is the dynamics of the growth. The evolution of the clus- 
ter ensemble is characterised by the number density Us of clusters of the mass s. 
For the simulated spectra of the Brownian coagulation we observed a self-similar 
evolution. Foii^juifh growth processes the mass spectrum can be described by the 
scaling ansat?^" 



[tf 



n,(i)^t— s-V(sA^), 



(7) 



where w, r, and z are scaling exponents and f{x) is the scaling function of the 
growth process. In the asymptotic limit, the scaling exponents are related by 
w = (2 — t)z. Furthermore, the mean cluster mass S{t) ~ YlT ''^sS^ evolves 
asymptotically as S{t) ~ t^ while the cluster nuniher density n{t) ~ ^^ Ug scales 
as n{t) r^ t~^ for r < 1 and n{t) ^ t~'^ otherwisecj. /^From the averaged spectra 
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Figure 6. Normalised evolution of the mean cluster size and the cluster number density for initial 
monomer number densities no of 2-10^m~^ (full line), 2-10*m~^ (dotted line), and 2-10^m~'^ 
(dashed line). The dashed-dotted line cpi*responds to the predictions based on the coagulation 
theory for the ballistic growth of fractals |12I. tii = {Kiino}~^ is the mean collision time between 
monomers. 



we calculated the mean cluster size and the cluster number density for different 
initial number densities no (Fig- h)- The calculated scaling exponents are related 
to each other as predicted by this concept. Consequently, the scaling ansatz ^ is 
appropriate to describe the dust growth driven by Brownian nrotion at low number 
densities. The scaling exponents z derived from the asymptotic growth of the 
mean cluster mass for different no are smaller than 2. Consequently, the fastest 
increase of the mean friction time due to Brownian growth is (tj) ~ t^'^. This 
result supports the point made in the previous section that the short time scales 
of dust growth required by astronomical and geophysical constraints can not be 
due to coagulation driven by Brownian motion. 

In numerical studies of the cosmic dust growth the aggregation process is 
very often described by a rate equation approached. As the spatial dependence 
of the number densities is neglected, this approach is a mean field description of 
the growth. In this framework the scaling exponents z, w, and r are ruled by 
the scaling behaviour of the rate coefficients and do not depend on the initial 
number density no. However, as indicated in Fig. o, our simulation results show 
a significant dependence upon no, which might be due to spatial fluctuations of 
the number densities. Our results suggest that using the Smoluchowski theory 
should be done with caution for describing the dust growth in proto-planetary discs. 



4 Conclusions 

In the early stages of dust growth the coagulation of dust particles driven by 
Brownian motion plays an important role. In order to explain the short time 
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scales of dust growth as constrained by astronomical and geophysical studies not 
only the mass of the dust particles but also their friction time has to increase 
quickly. We developed an N-particle model that allowed us to study coagulation 
driven by Brownian motion. In contrast to earlier studies, the dust growth 
was investigated at astrophysically relevant small number densities. Under such 
conditions the mean particle distance exceeds the typical aggregate diameter by 
orders of magnitude, and a collision will be a rare event. We derived a criterion 
which allows an efficient detection of candidates for imminent collisions. The 
N-particle-method is based upon an adaptive time step scheme respecting the 
individual dynamical states of the aggregates. Its basic concept is to perform 
on average constant "length steps", instead of using constant time steps. In 
order to minimise the influence of the decreasing number of particles within the 
simulation box, a new rescaling method is used throughout the aggregation process. 

We found that the structural properties of a dust cloud is characterised more 
completely by the distribution of the fractal dimensions of its constituent clusters. 
The mean of this distribution of D^ ~ 1.8 is in good agreement with the results 
of Monte-Carlo simulations of the ballistic Cluster-Cluster Aggregation (BCCA). 
The friction time of the aggregates growing due to Brownian coagulation scales 
with the mass as (t/) ^ m^/^°. 

In addition, the dynamical evolution of the dust ensemble has been investigated. 
We found that the mass spectrum of the ensemble evolves self-similar. For that rea- 
son, the mass spectrum of the dust growth due to Brownian motion at low number 
densities can be described by a scaling ansatz. We derived from the simulated 
spectra that the mean friction time of the dust ensembles increases with time as 
(tj) ^ t^/^. This indicates that the fast dust growth in proto-planetary is not due 
to the Brownian coagulation. 
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